tfd_f(ind) = c_d0_f(ind)+g_d0_f(ind)+exports0_f(ind);
tfs_f(ind) = l0_f(ind)+sum(ind2,iof0_f(ind2,ind));

parameter
aa_f(ind,ind2)
;

aa_f(ind,ind2) = iod0_f(ind,ind2);

mm = 0;
colflg = 0;
rowflg = 0;
mflg = 0;
iflg = 0;


while( (mm<500) and (mflg < 1),
mm = mm+1;

* COLUMNS
iflg = 0;

chk_c(ind2) = tfd_f(ind2) + sum(ind,aa_f(ind2,ind));
sum_c(ind2) = tfs_f(ind2) + sum(ind,aa_f(ind,ind2));

scal_c(ind2) = sqrt(chk_c(ind2)/sum_c(ind2));

loop(ind2,
if (  (abs(sum_c(ind2)-chk_c(ind2)) >= 0.00000000001),
iflg = iflg + 1;
);
);

loop(ind,
aa_f(ind,ind2) = scal_c(ind2)*aa_f(ind,ind2);
);

if( (iflg = 0),
colflg = 1;
);

iflg = 0;
chk_r(ind) = tfs_f(ind) + sum(ind2,aa_f(ind2,ind));
sum_r(ind) = tfd_f(ind) + sum(ind2,aa_f(ind,ind2));
scal_r(ind) = sqrt( chk_r(ind)/sum_r(ind) );

loop(ind,
if (  (abs(sum_r(ind)-chk_r(ind)) >= 0.0000001),
iflg = iflg + 1;
);
);

loop(ind2,
aa_f(ind,ind2) = scal_r(ind)*aa_f(ind,ind2);
);

if( (iflg = 0),
rowflg = 1;
);

* Check for convergence

if( (rowflg = 1) and (colflg = 1),
mflg = 1;
);


);


display
iflg
scal_c
scal_r
colflg
rowflg
mm
mflg
aa_f
;

ex_f(ind) = (tfd_f(ind) +  sum(ind2,aa_f(ind,ind2)))  - (tfs_f(ind)+sum(ind2,aa_f(ind2,ind)));

iod0_f(ind,ind2) = aa_f(ind,ind2);

